Do you believe your error bars?

Computing uncertainty estimates from Monte Carlo data usually assumes the data is normal (or can be sufficiently reblocked to be normal, according to the central limit theorem). But what if we don't have quite that much data - what effect would the underlying distribution have on the uncertainty estimates (error bars)?

Mervlyn Moodley looks at this question in the paper The Lognormal Distribution and Quantum Monte Carlo Data.

The interesting figures seem to be figures 3 and 9. In Figure 3, we see the confidence intervals change as the data is reblocked (and hence comes closer to normal). Figure 9 is generated from a real simulation, and it seems there is very little difference between the error bars computed assuming normality and those taking the underlying lognormal distribution into account.

Bilinear QMC

I'm reading a recent paper by Arias de Saavedra and Kalos titled "Bilinear diffusion quantum Monte Carlo methods" (PRE 67, 026708). They present an algorithm that uses a pair of walkers (the bilinear part) to sample the square of the ground state wavefunction (rather than first power of the wavefunction like DMC).

They claim it's useful for computing unbiased expectation values (see my post from Monday, Feb. 24 - although I'm not sure it would help with derivative operators) and energy differences.

The part I find most intriguing is the way they use importance sampling - they get reasonable results for hydrogen without it. And when they do use importance sampling, it's only to remove singularities due to cusps. I would guess that any successful scheme (ie, able to scale reasonably with system size) needs to use as an accurate a trial wavefunction as possible to guide the sampling.

The more I study the paper, the less I understand it. I guess the first step is to try reproduce their results for the 1-D harmonic oscillator.

Steal this research idea

Automatic control of Monte Carlo simulations - I want to know how to do it, but I don't necessarily want to figure it out myself. :)

• Determining equilibrium and phase transitions.
  Typically the particles are started in a lattice or random configuration and let run a whlie to equilibrate. I want automatic detection of equilibrium, which is especially tricky near phase boundaries. Also near phase transitions, the system can suddenly change phases.
  
• Adjustment of parameters for optimum efficiency
  The key quantity affecting the efficiency is the correlation time, but it's very noisy. Is there some other quantity that we could use that would behave similarly, but with less noise? Since this will only affect the efficiency, approximations are quite okay. Like assuming some analytic form for the autocorrelation function with a small number of parameters and fitting to it?
  

  
  Also note that SGA-like algorithms can be used for control as well as optimization. The decay parameter reaches a constant to follow the system, rather than continuing to decrease, as in optimization.
  

• Better control of DMC population
  
  Occasionally, I have trouble with the number of walkers increasing rapidly. In my codes, this means it exceeds some maximum value and stops the program. Sometimes the number of walkers jumps to a large but stable value, and the energy becomes far too low and clearly wrong. The hard part is that population control introduces a bias, and can't be too intrusive.
  

  
  For these items, are there any concepts from control theory or signal processing that would be useful?

Dynamic Monte Carlo

Checking the arxiv.org preprint server, there was an article titled An Introduction To Monte Carlo Simulations Of Surface Reactions by A.P.J. Jansen, which describes Dynamic Monte Carlo methods (Kinetic Monte Carlo is similar). One important input is transition probabilities. On page 21, "Estimates of the error made using DFT for such systems are at least about 10 kJ/mol."

Can QMC do better?

VMC Optimization: SGA vs. fixed sample reweighting

In VMC, we parameterize the wave function and vary the parameters until the minimum energy (or variance) is reached. Fixed sample reweighting generates a number of configurations from the MC process, and then uses a least squares minimization algorithm and a reweighting estimate for the energy to minimize the variance of the energies.

The Stochastic Gradient Approximation (SGA) takes a noisy estimate of the gradient and moves in that direction ("downhill"). The amount it moves is scaled by a variable that decreases each iteration (usually 1/n). It is guaranteed to converge, but it may take a while.

I don't like the fixed sample reweighting (FSR) method on aesthetic grounds. SGA feels more in line with the "Monte Carlo spirit". FSR feels like it is taking very coarse steps: sample configurations for a while, stop and do a minimization, repeat these two steps for a few iterations. Each step is large, and results in significant improvement. The SGA, on the other hand takes a large number of small steps. Each individual step doesn't give much improvement (or any improvement at all), but the steps accumulate to give a good answer.

Aesthetic considerations aside, FSR requires no analytic gradients for efficient implementation. SGA needs analytic gradients for most of the parameters (you could do numerical gradients, but the time gets much worse)

FSR does variance minimization, and SGA does energy minimization - there may be some differences in the resulting quality of the wave functions (JCP 112, 4935)

Derivative operators and forward walking in DMC

Using Diffusion Monte Carlo (DMC) to compute estimators other than the energy leads to a "mixed estimator"- an average that combines the trial wavefunction and the exact ground state (actually the fixed node ground state) - not what we really want. However, the number of children of each walker after a long time is proportional to the ratio of the fixed node ground state to the trial wavefunction. This leads to the forward walking procedure for obtaining correct averages over the fixed node ground state. See Casulleras and Boronat ( PRB 52, 3654) for examples, references, and implementation details.

Numerical values for derivative operators, such as the kinetic energy, are computed by applying the operator to the trial wavefunction. The forward walking procedure only corrects the distribution in the average, not the values of the estimator. (we really want the derivative applied to the exact wavefunction, not the trial wavefunction).

Computing the kinetic energy (and the pressure) by naively doing forward walking on the kinetic energy is incorrect. For the kinetic energy, we can use the total energy and the forward walked potential energy to obtain the kinetic energy.

Use Bennett's method for free energy calculations

That's the conclusion that Nandou Lu, Jayant Singh, and David Kofke reach in a recent paper (
JCP 118, 2977)

Bennett's method requires an estimate of the free energy, but that can be obtained via iteration. The authors argue that even without fully optimizing the solution, the free energy produced is still better than other methods based on overlap sampling.